package com.otp.webapp.action;

import java.io.Serializable;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

import com.otp.model.test.instance.TestInstance;
import com.otp.model.test.template.Item;
import com.otp.model.test.template.ResponseScoreOption;
import com.otp.model.test.template.TestTemplate;
import com.otp.service.util.MathematicalFunctions;
import com.otp.webapp.comparators.CompOptionsByObjectId;


public class AdaptiveMethods implements Serializable{

	private static final long serialVersionUID = 1L;
	private TestTemplate testTemplate;
	private TestInstance testInstance;
	private List<AnsweredItem> ansdItems;
	private Double abilityEAP=0.;
	private Double abilityML=0.;
	private Double abilityMAP=0.;
	private Double errorEAP;
	private Double errorML;
	private Double errorMAP;

	private static final int NQ=61;
	private static final double thetaMin = -3.;
	private static final double HPM=0;
	private static final double HPVAR=1.0;
	private static final int maxIterations=10000;
	private static final double NewtRaphError = 0.00001;


	public AdaptiveMethods(){
		abilityEAP=0.0;
		abilityML=0.0;
		abilityMAP=0.0;
		errorEAP=1000.0;
		errorML=1000.0;
		errorMAP=1000.0;
		ansdItems = new ArrayList<AnsweredItem>();
	}

	public Double getAbilityEst(){
		return abilityEAP;
	}	
	public Double getError(){
		return errorEAP;
	}
	public TestTemplate getTestTemplate() {
		return testTemplate;
	}
	public void setTestTemplate(TestTemplate testTemplate) {
		this.testTemplate = testTemplate;
	}
	public TestInstance getTestInstance() {
		return testInstance;
	}
	public void setTestInstance(TestInstance testInstance) {
		this.testInstance = testInstance;
	}
	public Double getAbilityEAP() {
		return abilityEAP;
	}
	public void setAbilityEAP(Double abilityEAP) {
		this.abilityEAP = abilityEAP;
	}
	public Double getAbilityML() {
		return abilityML;
	}
	public void setAbilityML(Double abilityML) {
		this.abilityML = abilityML;
	}
	public Double getAbilityMAP() {
		return abilityMAP;
	}
	public void setAbilityMAP(Double abilityMAP) {
		this.abilityMAP = abilityMAP;
	}
	public void setErrorEAP(Double errorEAP) {
		this.errorEAP = errorEAP;
	}
	public Double getErrorEAP() {
		return errorEAP;
	}	
	public void setErrorML(Double errorML) {
		this.errorML = errorML;
	}
	public Double getErrorML() {
		return errorML;
	}
	public void setErrorMAP(Double errorMAP) {
		this.errorMAP = errorMAP;
	}
	public Double getErrorMAP() {
		return errorMAP;
	}	
	public List<AnsweredItem> getAnsdItems() {
		return ansdItems;
	}

	public void setAnsdItems(List<AnsweredItem> ansdItems) {
		this.ansdItems = ansdItems;
	}

	public void addItem(AnsweredItem itm){
		ansdItems.add(itm);
	}

	public void updateDichot(){
		EAP(0);
		DichotomousMLandMAP();
	}

	public void updatePolyt() throws Exception{
		EAP(1);
		PolytomousMLandMAP();
	}



	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//DICHOTOMOUS MLE AND MAP USING NEWTON-RAPHSON														//
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	private void DichotomousMLandMAP(){
		Double f_thetaML=0.;//Equation to be solved f_theta=d(lnL)/dtheta=0
		Double fprm_thetaML=0.;//Derivative of f_theta, required by the Newton-Raphson method
		Double f_thetaMAP=0.;
		Double fprm_thetaMAP=0.;

		int iter=0;

		while(true){
			f_thetaML=0.0;
			fprm_thetaML=0.0;
			f_thetaMAP=0.0;
			fprm_thetaMAP=0.0;
			for(int i=0; i< ansdItems.size(); i++){
				Item itm = ansdItems.get(i).getItem();
				Double score = ansdItems.get(i).getScore();//Score must be either 0 or 1 !!!!
				Double a = (Double) itm.getA();
				Double b = (Double) itm.getB();
				Double c = (Double) itm.getC();

				//Index 1 corresponds to ML
				Double phat1=1.0/(1.0+java.lang.Math.exp(-a*(getAbilityML()-b)));
				if(phat1<0.00000001) phat1=0.00000001;
				if(phat1>0.99999999) phat1=0.99999999;
				Double pt1 = c+(1.0-c)*phat1;
				if(pt1<0.00000001) pt1=0.00000001;
				if(pt1>0.99999999) pt1=0.99999999;
				Double psp1 = phat1/pt1;
				f_thetaML+=a*(score-pt1)*psp1;
				fprm_thetaML+=a*a*score*psp1*(1.-psp1)-a*a*phat1*(1.-phat1);

				//Index 2 corresponds to MAP
				Double phat2=1.0/(1.0+java.lang.Math.exp(-a*(getAbilityMAP()-b)));
				if(phat2<0.00000001) phat2=0.00000001;
				if(phat2>0.99999999) phat2=0.99999999;
				Double pt2 = c+(1.0-c)*phat2;
				if(pt2<0.00000001) pt2=0.00000001;
				if(pt2>0.99999999) pt2=0.99999999;
				Double psp2 = phat2/pt2;
				f_thetaMAP+=a*(score-pt2)*psp2;
				fprm_thetaMAP+=a*a*score*psp2*(1.-psp2)-a*a*phat2*(1.-phat2);

			}
			f_thetaMAP-=(getAbilityMAP()-HPM)/HPVAR;
			fprm_thetaMAP-=1.0/HPVAR;

			Double delta1 = -(f_thetaML)/fprm_thetaML;
			if(delta1>0.5) delta1=0.5;
			if(delta1<-0.5) delta1=-0.5;
			abilityML+=delta1;
			if(abilityML<thetaMin) abilityML=thetaMin;
			if(abilityML>-thetaMin) abilityML=-thetaMin;

			Double delta2 = -(f_thetaMAP)/fprm_thetaMAP;
			if(delta2>0.5) delta2=0.5;
			if(delta2<-0.5) delta2=-0.5;
			abilityMAP+=delta2;
			if(abilityMAP<thetaMin) abilityMAP=thetaMin;
			if(abilityMAP>-thetaMin) abilityMAP=-thetaMin;


			if(iter++>maxIterations) {System.out.println("Max iterations in Newton Raphson reached !!!"); break;}
			if((java.lang.Math.abs(delta1)<NewtRaphError)&&(java.lang.Math.abs(delta2)<NewtRaphError)) break;
		}
		Double info1=0.0;
		Double info2=0.0;
		for(AnsweredItem ans:ansdItems){
			info1+=MathematicalFunctions.FisherInfo3PL(abilityML, ans.getItem().getA(), ans.getItem().getB(), ans.getItem().getC());
			info2+=MathematicalFunctions.FisherInfo3PL(abilityMAP, ans.getItem().getA(), ans.getItem().getB(), ans.getItem().getC());
		}
		errorML= 1.0/java.lang.Math.sqrt(info1);
		errorMAP= 1.0/java.lang.Math.sqrt(info2);
	}


	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
	//			POLYTOMOUS MLE AND MAP USING NEWTON-RAPHSON														//
	/////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
	private void PolytomousMLandMAP() throws Exception{
		Double f_thetaML=0.;//Equation to be solved f_theta=d(lnL)/dtheta=0
		Double fprm_thetaML=0.;//Derivative of f_theta, required by the Newton-Raphson method
		Double f_thetaMAP=0.;
		Double fprm_thetaMAP=0.;

		int iter=0;

		while(true){
			f_thetaML=0.0;
			fprm_thetaML=0.0;
			f_thetaMAP=0.0;
			fprm_thetaMAP=0.0;

			for(AnsweredItem ans:ansdItems){
				List<ResponseScoreOption> ansOptions = ans.getAvailableOptionEntities();
				Collections.sort(ansOptions, new CompOptionsByObjectId());//This is important
				Double a=ans.getItem().getA();
				List<Double> b=ans.getThresholds();
				Integer answer=ans.getSelectedOptionIndex();

				//Index 1 is ML
				//Index 2 is MAP
				Double grm1 = 0.;
				Double grmPrm1 = 0.;
				Double grmDblPrm1 = 0.;
				Double grm2 = 0.;
				Double grmPrm2 = 0.;
				Double grmDblPrm2 = 0.;

				grm1=MathematicalFunctions.GRM(answer, abilityML, a, b);
				grmPrm1=MathematicalFunctions.GRMPrm(answer, abilityML, a, b);
				grmDblPrm1=MathematicalFunctions.GRMDblPrm(answer, abilityML, a, b);

				grm2=MathematicalFunctions.GRM(answer, abilityMAP, a, b);
				grmPrm2=MathematicalFunctions.GRMPrm(answer, abilityMAP, a, b);
				grmDblPrm2=MathematicalFunctions.GRMDblPrm(answer, abilityMAP, a, b);

				f_thetaML+=grmPrm1/grm1;
				fprm_thetaML+=(grm1*grmDblPrm1-grmPrm1*grmPrm1)/(grm1*grm1);

				f_thetaMAP+=grmPrm2/grm2;
				fprm_thetaMAP+=(grm2*grmDblPrm2-grmPrm2*grmPrm2)/(grm2*grm2);
			}

			f_thetaMAP-=(abilityMAP-HPM)/HPVAR;
			fprm_thetaMAP-=1.0/HPVAR;

			Double delta1 = -(f_thetaML)/fprm_thetaML;
			if(delta1>0.5) delta1=0.5;
			if(delta1<-0.5) delta1=-0.5;
			abilityML+=delta1;
			if(abilityML<thetaMin) abilityML=thetaMin;
			if(abilityML>-thetaMin) abilityML=-thetaMin;

			Double delta2 = -(f_thetaMAP)/fprm_thetaMAP;
			if(delta2>0.5) delta2=0.5;
			if(delta2<-0.5) delta2=-0.5;
			abilityMAP+=delta2;
			if(abilityMAP<thetaMin) abilityMAP=thetaMin;
			if(abilityMAP>-thetaMin) abilityMAP=-thetaMin;


			if(iter++>maxIterations) {System.out.println("Max iterations in Newton Raphson reached !!!"); break;}
			if((java.lang.Math.abs(delta1)<NewtRaphError)&&(java.lang.Math.abs(delta2)<NewtRaphError)) break;
		}

		Double info1=0.0;
		Double info2=0.0;
		for(AnsweredItem ans:ansdItems){
			List<ResponseScoreOption> ansOptions = ans.getAvailableOptionEntities();
			Collections.sort(ansOptions, new CompOptionsByObjectId());//This is important
			Double a=ans.getItem().getA();
			List<Double> b = ans.getThresholds();

			info1+=MathematicalFunctions.FisherInfoGRM(abilityML, a, b);
			info2+=MathematicalFunctions.FisherInfoGRM(abilityMAP, a, b);
		}

		errorML= 1.0/java.lang.Math.sqrt(info1);
		errorMAP= 1.0/java.lang.Math.sqrt(info2);
	}


	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
	//			EAP ESTIMATE																						//
	//			model 0: Dichotomous 3Pl																			//
	//			model 1: Plolytomous GRM																			//
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
	private void EAP(int model){

		Double[] theta = new Double[NQ];
		Double[] prior = new Double[NQ];
		Double[] logLikelihood = new Double[NQ];
		Double[] posterior = new Double[NQ];

		theta[0]=thetaMin;
		prior[0]=MathematicalFunctions.normalPDF(theta[0], HPM, HPVAR);
		logLikelihood[0]=0.0;
		Double interval = -2*thetaMin/(NQ-1.);
		for(int k=1; k<NQ; k++){
			theta[k]=theta[k-1]+interval;
			prior[k]=MathematicalFunctions.normalPDF(theta[k], HPM, HPVAR);
			logLikelihood[k]=0.0;
		}

		for(int k=0; k<NQ; k++){
			for(AnsweredItem ans:ansdItems){
				if(model==0){
					Double a=ans.getItem().getA();
					Double b=ans.getItem().getB();
					Double c=ans.getItem().getC();
					Double score=ans.getScore();//Dichotomous scoring only!!!!
					Double p=c+(1.0-c)/(1.0+java.lang.Math.exp(-a*(theta[k]-b)));
					if(java.lang.Math.abs(score-1.)<0.000001) logLikelihood[k]+=java.lang.Math.log(p);
					else logLikelihood[k]+=java.lang.Math.log(1-p);
				}

				if(model==1){
					List<ResponseScoreOption> ansOptions = ans.getAvailableOptionEntities();
					Collections.sort(ansOptions, new CompOptionsByObjectId());
					Double a=ans.getItem().getA();
					List<Double> b = ans.getThresholds();
					Integer answer = ans.getSelectedOptionIndex();

					Double p = MathematicalFunctions.GRM(answer, theta[k], a, b);
					if(p==0.0) p=0.0000000001;
					if(p==1.0) p=0.9999999999;

					logLikelihood[k]+=java.lang.Math.log(p);
				}
			}
		}

		Double snum=0.0;
		Double snumSq = 0.0;
		Double sdem=0.0;

		for(int k=0; k<NQ; k++) {
			posterior[k]=java.lang.Math.exp(logLikelihood[k])*prior[k];
			snum+=theta[k]*posterior[k];
			snumSq+=theta[k]*theta[k]*posterior[k];
			sdem+=posterior[k];
		}

		//Debugging
		//		System.out.println("prior");
		//		for(int k=0; k<NQ; k++){
		//			System.out.println(theta[k]+"\t"+prior[k]);
		//		}
		//		
		//		System.out.println("Likelihood");
		//		for(int k=0; k<NQ; k++){
		//			System.out.println(theta[k]+"\t"+java.lang.Math.exp(logLikelihood[k]));
		//		}
		//		
		//		System.out.println("Posterior");
		//		for(int k=0; k<NQ; k++){
		//			System.out.println(theta[k]+"\t"+posterior[k]);
		//		}

		snum=snum/sdem; //Normalize
		snumSq=snumSq/sdem; //Normalize
		abilityEAP=snum;
		errorEAP=java.lang.Math.sqrt(snumSq-snum*snum);
	}
}



	
	



	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
	//			POLYTOMOUS MLE, MAP AND EAP ESTIMATE USING QUADRATURES												//
	//////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
	/*private void EAPPolytomousEst(){
		int NQ =61;
		Double a;
		int answer=-1;
		Double p;

		Double[] theta = new Double[NQ];
		Double[] ak = new Double[NQ];//This is the prior function
		Double[] lxk = new Double[NQ];//This the (not normalized) posterior distribution

		Double snum=0.0;
		Double sdem=0.0;
		Double la;
		Double xla;

		theta[0]=-3.;
		ak[0]=MathematicalFunctions.normalPDF(theta[0], 0.0, 1.0);
		lxk[0]=0.0;
		Double deltaTheta = 6.0/(NQ-1);
		for(int k=1; k<NQ; k++) {
			theta[k]=theta[k-1]+deltaTheta;
			ak[k]=MathematicalFunctions.normalPDF(theta[k],0.0,1.0);
			lxk[k]=0.0;
		}

		for(int k=0; k<NQ; k++){
			for(int i=0; i<ansdItems.size();i++){
				AnsweredItem ansObject = ansdItems.get(i);
				List<ResponseScoreOption> ansOptions = ansObject.getAvailableOptionEntities();
				a=ansObject.getItem().getA();
				List<Double> b = new ArrayList<Double>();
				for(int j=1; j<ansOptions.size(); j++) b.add(ansOptions.get(j).getDifficultyThreshold());

				//Collections.sort(ansOptions);//THIS IS VERY IMPORTANT!!!
				//

				for(int j=0;j<ansOptions.size(); j++){
					if(ansObject.getSelectedOptionEntity().getObjectId()==ansOptions.get(j).getObjectId()) {answer=j; break;}
				}
				p = MathematicalFunctions.GRM(answer, theta[k], a, b);
				if(p==0.0) p=0.0000000001;
				if(p==1.0) p=0.9999999999;

				lxk[k]+=java.lang.Math.log(p);
			}
		}

		Double MAPmax=-1000.;
		Double MLmax=-1000.;
		Integer MAPi=-1;
		Integer MLi=-1;
		for(int k=0; k<NQ; k++){
			if(lxk[k]>MLmax) {MLi=k; MLmax=lxk[k];} 
			lxk[k]=java.lang.Math.exp(lxk[k]);
			la=lxk[k]*ak[k];
			if(la>MAPmax) {MAPi=k; MAPmax=la;}
			xla=theta[k]*la;
			snum+=xla;
			sdem+=la;//Summing this gives the normalizing constant
		}
		abilityMAP = theta[MAPi];
		abilityML = theta[MLi];

		//This is the distribution
		for(int k=0; k<NQ; k++){
			Double y = lxk[k]*ak[k];
			//System.out.println(theta[k]+"\t"+y);
		}

		abilityEst=snum/sdem;
	}*/

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
//			VERY SIMPLE ESTIMATE																				//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////	

/*private void verySimpleEst(){

	double adminCount=0.5;
	double rightCount=0.5;
	double wrongCount=0.5;
	double sumDiffic=0.0;
	double sumSqDiffic=0.0;
	double meanDiffic=0.0;
	double disperseFactor=0.0;

	for(int i=0; i<ansdItems.size(); i++)
	{
		adminCount++;
		double score = (double) ansdItems.get(i).getScore();
		double maxScore = (double) ansdItems.get(i).getItem().getMaxScore();
		rightCount+=score/maxScore;
		wrongCount+=1-score/maxScore;
		sumDiffic+=ansdItems.get(i).getItem().getB();
		sumSqDiffic+=ansdItems.get(i).getItem().getB()* ansdItems.get(i).getItem().getB();
		meanDiffic = sumDiffic/adminCount;
		disperseFactor = java.lang.Math.sqrt(1+(sumSqDiffic/adminCount+meanDiffic*meanDiffic)/2.89);
	}
	abilityEst = sumDiffic+disperseFactor*java.lang.Math.log(rightCount/wrongCount);
	stdError = disperseFactor * java.lang.Math.sqrt(adminCount/(rightCount*wrongCount));
}*/